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Abstract. - In this study, the phase field model of crack propagation is used to study the dynamic 
branching instability in the case of inplane loading in two dimensions. Simulation results are in 
good agreement with theoretical predictions and experimental findings. Namely, the critical speed 
at which the instability starts is about 0.48cs . They also show that a full 3D approach is needed to 
fully understand the branching instability. The finite interface effects are found to be neglectable 
in the large system size limit even though they are stronger than the one expected from a simple 
one dimensional calculation. 
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Introduction. — The study of crack propagation has 
gained a lot of interest in the physical community dur- 
ihg the recent years [1]. Among the many reasons that 
have driven this interest is the difBculty to understand 
the mechanism leading to the branching instability [2, 3] 
that prevents a crack from reaching its theoretical limiting 
Speed: the Rayleigh wavespeed. During the branching in- 
stability, a single straight propagating crack separates into 
two sub cracks. While, the instability takes place in the 
process zone (the microscopic region ahead of the crack tip 
where the interatomic bond breaking leading to the crack 
propagation occurs), its effects (namely the propagation 
of two sub cracks) are macroscopic. 

, Hence a full description of the branching instability 
lieeds to describe both the macroscopic and the micro- 
scopic scales. The classical theory of crack propagation, 
describes well the macroscopic scale where the linear elas- 
ticity theory is valid but describes the crack propagation 
with laws that postulate the evolution of the crack path 
is a function of the stress intensity factors (number de- 
scribing the singularity at the crack tip, that correspond 
to the three modes of crack propagation (sec fig. 1). Such 
laws that do not describe in any way the small scale of the 
process zone can not account for the onset of the branch- 
ing instability unless one explicitly provides a criteria that 
determines when a crack divides into two sub-cracks. Yet, 
they provide useful results such as necessary conditions for 
branching [4,5]. 

On the opposite, the use of molecular dynamics simula- 



tions can help understanding the way a crack propagates 
and eventually divides into branches, since it allows to com- 
pute what happens at the microscopic scale. [6] Nonethe- 
less, it does not allow to reach neither long enough time 
scales nor large enough space scales to fully simulate mul- 
tiscale problems. 

Therefore, in order to describe crack propagation an al- 
ternative approach may be to use phenomenological mod- 
els of crack propagation. Such models aim at describing 
the behavior of the elastic material in the process zone 
with the use of non-linear elasticity. Typically, these mod- 
els describe the process zone as a region where the material 
softens gradually from a fully intact state where the law 
of elasticity are verified to a fully broken state where the 
crack cannot transmit any stress. 

A pioneering work was the introduction of the cohesive 
zone model where the crack line is prolongated by a soft- 
ening segment [7,8]. More recently using idea from the 
phase transition theory, phase-field models of crack prop- 
agation were introduced [9-16] to describe crack propa- 
gation. They have already proved they can well describe 
the complex crack path in many cases [10, 17]. In addition 
such models have been able to reproduce the feature of 
the crack branching in the case of mode III cracks [12, 18]. 
Here, a study of the branching instability using the phase- 
field model is presented in the case of the inplane crack 
propagation (mode I and II). It must be emphasized here 
that while previous works have been devoted to the case 
of mode HI crack propagation (paper tearing mode), this 
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Fig. 1: left The three difFerent modes of crack propagation: (a) 
mode III, (b) mode I and (c) mode II. rightCrack tip shape 
for difFerent crack speeds. Note the increase of the tip radius of 
curvature as the crack speed increases (• ii = 0.12, • v — 0.33 
and • 11 = 0.47 to compare with the branching threshold of 
0.52). For the liigher speeds, the crack surface is not exactly 
straight and presents some small irregularities. 



works is devoted to mode I crack propagation (pure ten- 
sile loading) which corresponds to usual experimental se- 
tups. The convergence of the numerical results is discussed 
and a comparison with available experimental data is pre- 
sented. Results show that the phase field approach is in 
good quantitative agreement with both theoretical predic- 
tions [5] and experimental results, they also show that in 
order to fully understand three dimensional crack branch- 
ing, three dimensional simulations are needed. 

model and numerical implementation. The 

model used here has been introduced in the case of an- 
tiplane loading in [9] It has been shown to reproduce well 
the branching instability of cracks under antiplane load- 
ing [18]. The extension of the model to crack propagation 
under inplane loading was presented in [10] where it was 
shown to reproduce well the oscillating instability of cracks 
under bi-axial strain and also some feature of the branch- 
ing instability. This later model was designed to comply 
with the principle of local symmetry and it has recently 
been used in the study of quasistatic crack propagation 
under thermal load. [16] 

The principle of the phase field approach is to introduce 
an additional continuously (in space and time) varying 
variable that will describe the state of the elastic mate- 
rial (e.g. (/) = 0: broken and = liintact). It will evolve 
obeying an evolution equation which couples (/> with the 
elastic fields. The equations of motion of the elastic mate- 
rial are also coupled in an appropriate way with the phase 
field. In this framework, the crack surface can be seen 
as an isocontour of the phase field (e.g. (/) = 0.5). Since 
(^ is varying continuously, the singularity at the crack tip 
is smeared out and the sharp tip is replaced by a smooth 
region where (j) is varying rapidly from 1 to 0. This re- 
gion may be considered as the process zone where failure 
occurs. This approach differs from the cohesive zone ap- 
proach [7, 8] because instead of smoothing the singularity 
by prolongating the crack line by a soft interface line di- 
rected along the predicted a priori path of the crack, it 
smooths the singularity at the crack tip by introducing a 



soft region. In this approach the crack faces are no longer 
lines but regions where the elastic moduli of the material 
go from their nominal values to zero and there is no need 
to know a priori the crack path. 

First the model is described in more details. The equa- 
tion for the displacement fields (ui where i is in x, y) 
writes as in [10]: 
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Eei — 



SE„ 
Sui 



(1) 



9m-itrie)r + nitrie')) (2) 



where e is the strain tensor and is, in the limit of linear 
elasticity, equal to eij = (diUj + djUi)/2. The function 
g = 4(/)^— S^** couples the phase field with the displacement 
fields and it has been chosen so that the stress transmitted 
through a crack vanishes in the large system limit (i.e. 
when the size of the system becomes much larger than 
the width of the crack interface introduced by the phase 
field.) [9]. The equation for the phase field writes: 
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Equation (3) insures that the phase field will not increase, 
which means that a crack can not heal. Equations (5, 4), 
are similar to the free energy used in [9] . The form of the 
coupling term in (eq. 4), indicates (if one neglects cur- 
vature effects at the tip) that the crack will propagate in 
a region where the value of f^ is higher than Cc- Con- 
trarily to the model of [9] , here £^ is not the local energy 
density. It takes into account the fact that the material 
is either under compression (ir(e) < 0, ^^ is the contri- 
bution of shear to the local elastic energy density) or ex- 
tension (ir(e) > 0, £^ is the local elastic energy density) 
and that a crack should not propagate in a material under 
compression. Therefore equations (5, 4) impose that in 
a material under compression, the favored state is intact 
material (0=1) 

This model depending on the imposed boundary condi- 
tions can describe the propagation of a crack (or of many 
cracks) under arbitrary loading conditions. Model param- 
eters used here are ec(= 1), D{= 0.25), /3(= 1, 2, 4 and 8) 
and (5(= 0.1). The fracture energy in the large system 
[9, 10] limit is: 



r = 




VsvW+W^gW) 



(6) 



Here we have applied this model to the study of a dy- 
namic crack propagating under mode I (tensile) loading: 
the top and bottom boundaries of the sheet are moved of 
Uy{±W/2) = ±Ay/2 and Ux{x = 0/L) = 0. Parameters 
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arc D = 0.25 so that the thickness of a crack (w^) at the 
onset is approximately 1 s.u.(I? = 1 was also used with- 
out any significant change). A = 1, /i = 1, which implies 
that the shear wave speed is Cs = 1 and the Rayleigh wave 
speed is Cr = 0.91. The system is prepared at mechanical 
equilibrium with an initial prc-crack. Then the equations 
(1,3) are simulated on a regular grid (with grid spacing 
dx ~ 0.1) using a forward Eulcr scheme. The discretiza- 
tion scheme is chosen to keep the mechanical energy con- 
stant (if the phase field is kept constant and so that the 
total energy decreases when the phase field evolves.). In 
addition, checks were performed to insure that discretiza- 
tion effects are ncglectable: dividing the grid spacing by 2 
did not affect results by more than 1 percent. 

Since the focus of this study is the stability/instability 
of straight cracks propagating at constant speed, simula- 
tions were performed on a grid moving with the crack tip 
in order to let the transient accelerating regime disappear. 
The simulation grid had an aspect ratio [L/W) of 1. Sim- 
ulations using larger and smaller aspect ratios (0.5, 2 and 
4) showed a small change in quantitative results: namely 
an approximately drop of the branching speed of 2 per- 
cents in the infinite aspect ratio limit for the three smallest 
width (12.5 25 and 50 and 100) and /3 = 1. On the con- 
trary, changing the ratio W/w^, where w^ is the interface 
thickness of the crack, did affect quantitatively the results. 
But, as expected [9], results converge towards a limit when 
the size of the sample W . compared to the width of the 
crack w^, goes to infinity (keeping the clastic energy stored 
in the micrackcd medium ((A/2 + /i)L(A.y/L)'^) per unit 
length and other parameter constant). Therefore the fact 
that results are affected when changing W does not indi- 
cate that the branching is size dependent. It is indicative 
of the effects of the finite width of the crack in the phase 
field model (It is not due to the fact that the full stress 
relief only occurs in the infinite W limit, which is a ID 
effect. It is due to the fact that a finite distance exists 
between the edges of the crack). 

Straight crack. — First some results on the prop- 
agation of a stably propagating single crack arc pre- 
sented. Using the model parameters, the fracture energy 
is 2r — 0.971. According to the Griffith criterion, one 
expects that a crack will begin to propagate when the 
elastic energy stored per unit length is larger than the 
energy needed to create two crack interfaces, that is for 
Ay = yJ{TW)/{\/2 -\- //). Simulations results show that 
this prediction is extremely well verified, with an error 
smaller than a percent, for the smallest simulation boxes 
used here (W=25). For bigger systems (W=50, 100 and 
200), the error is always at most of the same order of 
magnitude. Hence the model used here is in very good 
agreement with the Griffith criterion. 

When the value of Ay is increased, the speed of the 
straight crack v is expected to obey the law: 
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Fig. 2; leftCrack speeds as a function of the available elas- 
tic energy (oc A^/VK) normalized by the crack energy (derived 
from the model equations as in [9]) for /3 = 1 and different sys- 
tem sizes. One should note that the curves collapse extremely 
well on a master curve and that, for the sake of simplicity, only 
the regime of straight crack propagation is shown. right Frac- 
ture energy normalized by the fracture energy at zero velocity 
computed using eq.7 as a function of crack speed for different 
values of /3 (/3 =1 (□),2 (*),4 (X) 8 (+)). The lines correspond 
to the law eq. 8 and fit well the numerical datas. For higher 
values of v one can observe a deviation from the linear law that 
corresponds to the onset of the branching instability. 

where T{v) is the mechanical energy that is either diss- 
pated or converted into surface energy when the fracture 
advances of one space unit at speed v, G is the mechanical 
energy available far ahead from the crack tip and Ai{v) 
is a universal function that depends solely on the elastic 
properties of the material (for further details see ref. [19] 
sections 5.3 and 5.4). Simulation results using the phase 
field model and eq. 7 show that 



T{v)^T{0){l + apPv), 



(8) 



GAi{v) = T{v) 



(7) 



with a/3 « 0.5 a constant, for values of v below the branch- 
ing threshold, and indicate that for /3 — > 0, the fracture 
energy is indcpendant of the crack speed (zero dissipation 
limit). It must be noted here that for or a given model pa- 
rameter set, the relation between the steady crack speed 
V and the available mechanical energy ahead of the crack 
tip is independant of the system size (see fig. 2). 

At this point, before turning to the study of branch- 
ing instability, it is worth to summarize the results ob- 
tained here. First the phase field model, as already shown 
in [10] obeys the Griffith criterion. In addition, the re- 
lation between the model parameter (3 and the speed de- 
pendance of the fracture energy has been computed (eq.8). 
The later result together with theoretical prediction in [5] 
which show that the branchinsg speed Vc behaves like 
a+&(F(iJc)-F(0))/F(0) , allows to predict that the branch- 
ing speed in the phase field model must behave linearly 
with /? (taking into account the fact that in [5], only a 
necessary condition for branching was derived while here 
we focus on the onset of branching). 

We now turn to the main purpose of this work: the 
study of the branching instability. 

branching instability. — When the value of Aj, is 
further increased, the crack tip becomes unstable and the 
crack divides into two cracks. As shown in [18], this in- 
stability occurs through tip blunting. As the crack speed 
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Fig. 4: Speed of the crack as a function of time. One can 
see the intermittent behavior of the crack speed. The burst 
of crack speed oscillation correspond to branching events while 
the regions where the speed is varying slowly correspond to a 
situation where the crack has divided into two branches that 
propagate together for a while before one of them takes over. 
l3 = 8. 



Fig. 3: Typical crack paths during the branching instability 
(the crack is propagating from left to right). The top plot 
shows a single branching event. The bottom plot shows a fast 
crack with multiple macroscopic branches appearing. In the 
top plot, the middle of the elastic plate is at y = 25, and in 
the bottom it is at j/ = 50. Parameter values are: /3 = 1 and 
G/Gc = 2.5 



increases, the radius of curvature of the crack tip increases 
and when the speed is larger than a critical value, the 
crack tip splits into two branches that propagate (see fig. 
3, where one can see that the crack tip is still much smaller 
than the system width). One of those two branches can be 
faster than the other and through screening of stress even- 
tually prevents the propagation of the other branch. Then 
if strain is sufficient, the scenario repeats and one can see a 
succession of branching events giving birth to two cracks, 
the selection of one of them its acceleration and separa- 
tion into two cracks... This regime of sustained disordered 
branching events can only take place if strain is sufficient 
and corresponds to the macro-branching instability (see 
fig. 3). Its signature when considering the crack surface is 
the presence of successive branches of various sizes. In ad- 
dition to this qualitative description, one notices a sharp 
increase in the variance of the crack speed which is similar 
to the one observed in experiments despite some discrep- 
ancies (The variance of the crack speed as a function of 
the mean crack speed increases via a discontinuous jump 
in 2D simulations while in experiment it increases linearly 
above threshold) which could be attributed to the 3D na- 
ture of experiments. This transition is easily noticed when 
considering the velocity profile (Here crack velocity is the 
time derivative of the position of the most advanced point 
of the crack: the iso contour ~ 05). Instead of small 
oscillations one can see sharp decreases of the crack speed 
followed by a progressive acceleration of the crack and 
again a sharp decrease (see fig. 4) . This behavior of the 
crack speed is typical of a branching event and is in good 
agreement with the scenario presented in [5] , with numer- 
ical observations of [18] and with experimental results. 
The critical speed depends on the model parameters 



(here: f3^) and on the system size (see fig. 5). While the 
dependence on (3 is expected since /? measures the energy 
dissipation at the crack tip, the dependence on W/wtj, is 
not the result of any physical phenomena and comes from 
the way the phase field models deals with the singularity 
at the crack tip. 

Indeed, it does not allow to consider an infinitely sharp 
crack. It only allows to reproduce a crack of finite width 
Wtj, and one should retrieve the results of an infinitely sharp 
crack in the limit w^/W — *■ 0. Simulations results show 
that this is actually the case. In the limit w^/W -^ 0, 
the value of the critical speed at which the macroscopic 
branching instability occurs for a given value of f3 con- 
verges toward a finite limit while the qualitative behavior 
of the system remains unchanged. As mentioned earlier, 
this limit is a function of /3 which measures the dissipa- 
tion in the material. As a result, in the limit /3 ^ 0, the 
phase field model reproduces the behavior of a perfectly 
brittle material. Simulations using different values of /3 
show that the branching speed scales linearly with /3 and 
converges toward a finite value when /3 — > as one would 
expect from [5] and the results of the study of steady crack 
propagation. 

Then in the limit /3 ^ 0, {w^/W) -^ 0, one gets the 
value of a critical macroscopic branching speed for a per- 
fectly brittle material without any dissipation of approxi- 
mately 0.48. This value is in fairly good agreement (a 15% 
discrepancy) with theoretical results of [5] which predicts 
a branching speed of 0.42 (with parameters used here). 
The discrepancy can be easily attributed to the fact that 
in [5], the critical speed that was computed, was the speed 
at which branching can occur, i.e. a necessary condi- 
tion, while here the computed speed is the one at which a 
straight crack is no longer stable. And to reinforce this, it 
is important to mention that during simulations branching 
events were observed at instantaneous speeds lower than 
the the measured critical velocity. 



^Additional simulations have shown that dividing 5 by 10 did not 
affect results and simulations of [18] in the case of mode III cracks 
show that (5 — > 0. is the correct limit in the sense that it allows to 
retrieve theoretical predictions 
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Fig. 5: (a) critical speed for which the straight crack is no 
longer stable as a function of /? for different value of W/w^, 
the ratio of the system size (W) on the width of the phase field 
crack {w,/,. One should note that as /3 goes toward zero (zero 
dissipation limit), the critical speed converges toward a value 
that depends on the system size and that contrarily to what is 
observed in fig. 2 the curves do not collapse on a master curve. 
(x, +, triangles and circles correspond respectively to sizes of 
25, 50, 100 and 200). (b) Branching speed in the zero /9 limit 
(circles) and for /3 = 2 (diamonds) as a function of the system 
size W relative to the crack width (w,/,). The lines are used as a 
guide to the eye. The fact that the convergence of the speed is 
in \/'w^/W may be attributed to the effects of the finite width 
of the crack interface on the singular behavior at the crack tip. 
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Fig. 6: Typical curves of the crack speed as a function of time 
during repetitive branching events for two different values of /3 
(top: 8 and bottom: 4). One can see that the typical interval 
between maxima's of the (3 — 8 curve are approximately two 
times bigger than in the /3 = 4 curve. 



When comparing, the results obtained here with the 
experiments [20-22], one must first note that the three 
dimensional character of the branching instabihty cannot 
be reproduced by our simulations and that in PMMA, 
according to [22] (figs. 19 and 17), the branching insta- 
bihty can be described as 2 dimensional for average crack 
speeds much higher (w 0.58cj^) than the critical speed for 
the 3D branching instabihty (0.34cij). In addition, dis- 
sipative effects are present in PMMA and are not well 
quantified above the first critical speed (0.34c/j). Then, it 
is clear that our simulations can not be compared quanti- 
tatively with the experimental results. Nonetheless, some 
qualitative comparison are possible. First of all, the com- 
puted branching angle defined as the angle between the 
secondary crack and the main direction of propagation of 
the crack is in good agreement with experimental results. 
It is on average equal to 27° (with some visible variations), 
which is close to the experimental value of 30°. Interest- 
ingly enough, if one defines the branching angle as the 
angle between the two branches the computed values is 
then around 45° which is in rather good agreement with 
theoretical predictions of 54° [4,5]. The discrepancy in the 
computed values of angles may come from the fact that the 
phase field method does not allow to compute the branch- 
ing angle at the onset of the branches leading to a small 
under-estimation of its value. Indeed, the branching point 
is defined up to a few interface thicknesses and the behav- 
ior of the branches can only be described a few interface 
thicknesses away from the branching point. These short- 
comings of the phase field approach do not allow to check 
the asymptotic behavior of branches close to the onset of 
branching. Simulations in larger systems are needed to 
overcome the finite interface limitation. 

From a statistical point of view, the crack speed fluc- 



tuations power spectrum in the case where no dissipation 
occur in the bulk did not exhibit any clear frequency. The 
same behavior was observed in the case where a small 
dissipation was added in the bulk. Nonetheless, during 
the repetitive branching burst, the frequency of branch- 
ing events was found to be fairly constant for given model 
parameter values. Namely its Independence on the sys- 
tem size, on the crack speed could not be determined. 
Nonetheless, a clear dependence on the value of /? was 
found: the inter-branch interval is approximately propor- 
tional to /3 as can be observed in fig. 6. This behavior 
is true for small frequencies (i.e. for (3 > 2). But for 
small values of /3 (lower than 2), one can see that the 
inter-branch intervals ceases to decrease and saturate at 
a given value. This saturation may be attributed to the 
finite thickness of the interface (there exists a minimal in- 
terbranch distance which can be estimated to be equal to 
a few interfaces thicknesses. Here, assuming a crack speed 
of around O.bcs, and considering the crack width of 1. su., 
the highest frequency should be a fraction of 0.5.^). Hence 
as /3 is decreased, phase field simulations show an increase 
of the typical frequency until it reaches the saturation fre- 
quency. 

Discussion. — Results of simulations presented here 
show that the phase field approach is useful to understand 
fracture as a pattern formation mechanism. Nonetheless, 
the fact that crack propagation is due to the l/\/r singu- 
larity at its tip, finite size effects due to the phase field 
approach turn out to be much stronger than one would 
expect from the ID calculation of [9]. This slow conver- 
gence does not prevent the model to retrieve quantitatively 
good results in the large system size and the use of bet- 



^It cannot be more than 0.5 and it is reasonable to assume that 
one needs at least a few crack width between branches. 
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ter simulation methods (such as Finite elements or adap- 
tive remeshing) should allow to consider systems where 
finite size effects are neglectable. Hence, the phase field 
model is able the reproduce correctly the branching in- 
stability in the case of 2D crack propagation, predicting 
a critical speed of 0.48cs, which is 15% larger than the 
necessary condition derived in [5] . Since the branching in- 
stability is essentially a 3D process [22], comparison with 
experiments is only indicating that the phase field model 
can reproduce the main features of the branching insta- 
bility. But results and comparison with experiments in- 
dicate that full 3D simulations are needed to be able to 
understand the branching instability. One of the expected 
results of this is to be able to understand what are the 
parameter that govern the formation of fracture patterns 
such as the parabolic markings shown in [22]. 
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